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SUMMARY 

The estimation of potential fields such as the gravitational or magnetic potential at the 
surface of a spherical planet from noisy observations taken at an altitude over an in- 
complete portion of the globe is a classic example of an ill-posed inverse problem. 
Here we show that the geodetic estimation problem has deep-seated connections to 
Slepian's spatiospectral localization problem on the sphere, which amounts to find- 
ing bandlimited spherical functions whose energy is optimally concentrated in some 
closed portion of the unit sphere. This allows us to formulate an alternative solution to 
the traditional damped least-squares spherical harmonic approach in geodesy, whereby 
the source field is now expanded in a truncated Slepian function basis set. We discuss 
the relative performance of both methods with regard to standard statistical measures 
as bias, variance and mean-squai^e en^or, and pay special attention to the algorithmic 
efficiency of computing the Slepian functions on the region complementary to the axi- 
symmetric polar gap characteristic of satellite surveys. The ease, speed, and accuracy 
of this new method makes the use of spherical Slepian functions in earth and planetary 
geodesy practical. 

Key words: Geodesy, Satellite Geodesy, Spectral Analysis, Inverse Theory, Statistical 
Methods, Spherical Harmonics 



>1 INTRODUCTION 



Satellites mapping out the spatial variations of the gravita- 
tional or magnetic fields of the Earth or other planets ideally 
fly on polar orbits, uniformly covering the entire globe. Thus 
potential fields on the sphere are usually expressed in spher- 
ical harmonics, basis functions with global support. For var- 
ious, especially engineering, reasons, however, inclined or- 
bits are favorable. These leave a "polar gap": an antipodal 
pair of axisymmetric polar caps, typically less than 10° in 
diameter, without any data coverage. Estimation of spheri- 
cal harmonic field coefficients from an incompletely sam- 
pled sphere is prone to error, since the spherical harmonics 
are not orthogonal over the partial domain of the cut sphere. 

The historically somewhat neglected geodetic po- 
lar g ap problem has been re vived by, among oth- 
ers. ISneeuw & van GelderenI Jl997h . and recently, 
lAlberteUa et al.l ( Il999h . who constructed a new basis 
of so-called Slepian functions (after ISlepian 1983) on the 
sphere. These bandlimited functions are designed to have 
the majority of their energy optimally concentrated inside 
the latitudinal belt composed of the entire globe minus 
the polar gap, i.e. the region covered by satellites. Slepian 



functions are orthogonal on both the entire as well as the cut 
sphere, a property that can be exploited to our advantage. 
Here, we study the inverse problem of retrieving a potential 
field on the unit sphere from noisy and incomplete obser- 
vations made at an altitude above their source. We derive 
exact expressions for the estimation error due to the tradi- 
tional method of damped least-squares spherical harmonic 
analysis as well as that arising from a new approach using a 
truncated set of Slepian basis functions. 

We cast the geodetic estimation problem in the much 
wider context of spatiospectral localization, whereby band- 
limited functions are spatiallv concentrated to regions of 
arbitrarv shape on the sphere (IWieczorek & Simonsll2005t 
ISimons et al.l 120061) . and derive a new semi-analytical nu- 
merical method to calculate the spherical Slepian functions 
on the latitudinal belt or its complement, the double polar 
cap. Our approach requires no numerical integration, and 
avoids the construction of matrices other than a tridiago- 
nal matrix whose elements are prescribed analytically. Find- 
ing spherical harmonic expressions for bandlimited func- 
tions concenttated to polar caps or latitudinal belts, as in 
Figure n thus becomes so effortless as to be achievable by 
a handful of lines of computer code, and the problems with 
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numerical stabili ty known to plag ue alternative approaches 
JAlbertellaetal.lll999t iPail et all IMl are avoided alto- 
gether. 

The key to this "magic" lay hidden in two little- 
kn own studies published se veral decades ago: the work 
by iGilbert & SlepianI lll977h on doubly orthogonal poly- 
nomials, and that on c ommuting differential operators by 
iGriinbaum et alJ (Il982h . It must be remem bered that one 
of Slepian's main discoveries (see, e.g., ISlepianI Il983h 
was the existence of a second-order differential oper- 
ator that commutes with the spatiospectral localization 
kernel concentrating to intervals on the real line. Cast 
in matrix form, finding the prolate spheroidal functions 
amounts to the diagonalization of a simple tridiagonal ma- 
trix (see, e.g. Percival & Walden 1993). In their study, 
[Gilbert & Slepian ( 1977) presented two additional commut- 
ing differential operators, which are applicable to the con- 
centratio n of Legendre polynom ials to one-and two-sided 
domains. iGrtinbaum et alJ (ll982^ proved that the matrix ac- 
companying the locaUzation to the single polar cap is, once 
again, tridiagonal. Here, we show this is also the case for 
the antipodal double polar cap and its complement, the lati- 
tudinal belt. The tridiagonal matrix elements coding for the 
single pola r cap, and their sol utions, were published by us 
elsewhere dSimons et al 1 120061) . The expressions applicable 
to the double polar cap appear here for the first time. 

The problems posed and solved in this paper are not lim- 
ited to geodesy and observations made from a satellite. In 
geomagnetism, our observation level may be the Earth's sur- 
face, and the source level at or near the core-mantle bound- 
ary. In cosmology, the unit sphere constituting the sky is ob- 
served from the inside out, and the galactic plane masking 
spacecraft measurements has the shap e of a latitudinal belt 
llTegmarkllT996t iHinshaw et al.l l2003l). Ground-based astro- 
nomical measurements may be confined t o a sm all circular 
patch of the sky ( Peeble s 1973: Tegmark || 1995l) . Finally, in 
planetary science, knowledge of the estimation statistics of 
properties observed over mere portions of the planetary sur- 
face is important in the absence of groundtruthing observa- 
tions. 




Figure 1. Geometry of the geodetic estimation problem, and some symbols 
used in this paper. In the lower left, an axisymmetric polar cap, shaded, of 
colatitudinal radius 0. In the lower right, an antipodal pair of polar caps, 
shaded, representing the geodetic polar gap. 

ureQ illustrates the case in which the region is a double 
polar cap symmetric about the polar axis z. The angular ra- 
dius of the polar caps is denoted by O. The double polar cap 
is representative of the geodetic case in which R is the so- 
called polar gap of missing observations; its complement R 
is a latitudinal belt of angular width tt — 20 around the equa- 
tor, as shown. In the following, for brevity, we will shorten 
all double summations to a notation requiring only a single 
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2 STATEMENT OF THE PROBLEM 

We are concerned with estimating source-level potential 
fields from noise-contaminated satellite observations at an 
altitude over an incomplete portion of the unit sphere. The 
geometry of this problem is illustrated in Figure^ The unit 
sphere f2 on which the unknown signal is defined is parame- 
terized as usual in terms of spherical coordinates, colatitude 
9 and longitude 0. The angular distance between two posi- 
tion coordinates f = [6, (j>) and f' = {9' , (p') is denoted by 
A. In the lower right, the domain over which satellite ob- 
servations are available is left unshaded, whereas the area in 
which measurements are missing is shaded grey. We denote 
the white region covered by satellite tracks by R, and the 
shaded, uncovered region by R. Although our treatment will 
start out quite general, without restrictions on the shape of R 
or R, as long as they are complementary closed regions on 
the surface of the unit sphere, the lower right panel of Fig- 



2.1 Preliminary considerations on tlie source signal 

We model the geophysical signal as a broadband, square- 
integrable, real-valued function s(f) on the surface of the 
unit sphere il = (9, (p), defined by the transform pair 

s(f) = s;,„Yi„i(f), Sim = / s{r)Yimir) dn. (1) 

Im 

The integers I and m are the degree and order of the real 
spherical harmonics y;,„(f )■ These are defined by 

{V2Xim{9) COS if -/ < TO < 

Xm{9) if TO = (2) 

V2Xi,n{9)sinm(f> if < m < 1, 
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2-1)' 



(3) 



(4) 



where PimiiA is the associated Le gendre function 
alization 
10.1998.) 

{l-m)\ 



e gendre fu r iction , 

and the nor malization constant (e.g. lEdmondsl Il996t 
lDahlen&Tromal998.) 



Ci^ = (2^ + 1) 



(5) 



{l + m)V 

With these choices the harmonics Yi,n{0,(f>) are orthonor 
maHzed on the unit sphere: 



Ylm {r)Yi>m> {r)dQ^ 6w Sm-m' ■ (6) 



The fixed-order orthogonahty relation for Xim (0) is 

Xi,nXi>m sin ed0 ^ 5w ■ (7) 

ZTT 

The addition theorem expresses the sum over all orders of 
spherical harmonics at different positions in terms of the an- 
gular distance A = arccos(f • f ') between them as 



m — — l 



An 



(8) 



where we note that Piq = Pq and Pi{Q) = 1. The delta 
function 6{r, f') = {sme)-'^6{e - e')S{(t) - (j)') defined by 



1=0 



Aw 



has the usual sifting property 

6{r,r')f{r)dn^ f{r'). 



(9) 



(10) 



The sum over all degrees to infinity of the fixed-order colat- 
itudinal functions at different arguments results in the colat- 
itudinal delta function: 



{sme)-'6{e - e') = 2^ ^ xunmie'). 



(11) 



where the upward continued signal c oefficients are given in 
terms of the source-level terms by dStacevI (19921 iBlakelvl 
(19951) 



a 



-i-i 



Sir. 



(14) 



The data over the region of coverage R would be given by 
eq. ( I13t were it not that they are contaminated by noise. In 
the uncovered areas R, no measurements are available. A 
satellite thus observes 



d(f) 



S|-(f) + n(f) if f e i? 
unknown if f e ^. 



(15) 



We will restrict attention to the case in which the measure- 
ment noise 71(f) is additive and given by a zero-mean sto- 
chastic process, which we assume to be white: 



(n(f)) 
(n(f )n(f')) 



0, 

N6{r,r'). 



(16) 
(17) 



Thus, the power of the noise is denoted by TV, and we use 
angular brackets to denote the ensemble averaging over all 
possible realizations required to define the process mean and 
its spatial (co)variance. 

Combining eqs ( I13> -( I14> and Q with the definition ((Sj, 
we can write the signal observed at orbital level as a convo- 
lution of the surface-level signal in the form 



st(f)- / r(f,f').s(f') dn\ 

Jn 

where we have defined a "point spread function" 



r(f,f') = ^(l + a)-'-i ('^^)P,(f.f'). 

1=0 



An 



(18) 



(19) 



Thus, the value of the potential field that is observed at a 
point f outside the unit sphere is a weighted mixture of the 
function values at f and distant other points f ' on the unit 
sphere. Measurements taken by a satellite at a > are af- 
fected by regions it does not fly over directly. A satellite 
thus does probe into the uncovered regions; conversely, in 
regions of coverage, it may be affected by uncovered areas. 
As the fractional altitude a increases, the convolution ker- 
nel r(f , f ') is increasingly supported globally. On the other 
hand, when a ~ 0, eq. ( I19> returns the delta function, eq. (|9}, 
and eq. ilEl merely illustrates its sifting property dlOt . 



2.2 Noisy measurements at satellite altitude 

For convenience we separate the signal into a bandlimited 
portion restricted to the degrees I ^ ^ L and a portion 
over the degrees I = L + 1 ^ 00 that complement it: 

L 00 
S{r) ^^ShnYijnir) + ^ SijnYlm{r), (12) 

I'm lrn>L 

where we define L to be the spherical harmonic bandwidth. 
At the satellite altitude a above the unit sphere the analytic 
signal is given by 

L 00 
S|(f) = ^St/,„yi,„(f) + ^ SijlmYhn{r), (13) 
I'm lm>L 



2.3 A new basis for bandlimited field estimators 

We seek an estimate s(f ) of the signal in eq. Q, at the level 
of the source, from the data d{r), given by eq. ( !15> , at alti- 
tude a. It is crucial to realize that, although any real physical 
signal s(f ) will in general be infinite-band, our estimate s(f ) 
must always be bandlimited. We are thus at liberty to define a 
new, bandlimited set of basis functions, to replace the spher- 
ical harmonics. In this manner, the broadband source field 
can be expressed as 

{L + lf 00 

s{r) = ^ Saga{r) + ^ si„iYim{r). (20) 

a=l lm>L 
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whereas the bandhmited estimated field is given by 

L (L+lf 
S(f) = ^S;„iy;m(f) = ^ SQ5a(f)- 



(21) 



Ira 



a = l 



These new bandlimited basis functions (?Q(f) indexed by 
a ~ 1, . . . , (i + 1)^ will themselves be combinations of 
spherical harmonics, inasmuch as they are defined by the 
transform pair 

L 

ga{v) = ^ga.lrnXlrniY), (22a) 
Im 

{L+lf 

Yirn{v) = ^ ga.hn9a{r), (22b) 



ga,im = I ga{r)Yijnir)dn. (22c) 

The new basis will be rendered orthonormal by requiring 
that 



gair)g[3{r)dQ = Sap, 



ga,hnga,l'm' = S„ 



(23a) 
(23b) 

(23c) 



The transformation of the spherical harmonic basis coeffi- 
cients Sim of the estimate to the expansion coefficients Sa in 
the new basis is achieved by 

(L + lf L 

Sim = ^ 9a,lmSa and Sa ga,lm.Slm, (24) 

OL — l Im 

as can be easily deduced by combining eq. Mil with the or- 
thonormaUty conditions, eqs and j23> . and using eq. J22> . 

Upward continued to the satellite altitude a, the estimate 
in either basis is 



{L + lY 



(25) 



where the spherical harmonic coefficients are naturally given 
by 

S]lm = (1 + a)"'"^s;„, (26) 

and the upward continued expansion coefficients of the new 
basis by 



St 



Q = ^.g«,/7n(l + a) ' ^ 



Sir. 



(27a) 



lin 



= X] W^9a,ini{'^ + a)~'^~^gp.im\si3, (27b) 

/3=1 \ hn I 

as is verified by combining eqs J24l i and j26t . As expected, 
eq. J27l i reduces to a trivial identity at the surface of the unit 
sphere, i.e. when a = 0, by virtue of eqs ( I23t or ( I24> . 



Given that the measurements made by the satellite are re- 
stricted to the domain R on the unit sphere fi, we consider it 
natural to require of the new basis functions that they be op- 
timally concentrated on this domain. We will seek a basis of 
functions g(f ) whose energy is maximally concentrated in- 
side of the domain R by maximizing the spatial energy ratio. 



g^ dn 



= maximum. 



(28) 



g^ dn 



Eq. ( I28t is a statement of Slepian's pro blem, a clas- 
sic in one-di mensional ti me-series analysis ( ISlepianlll983L 
|Percival & Wa ldeniri993h . on the two-dimen sional sphere, 
which we have recently studied in detail (ISimons et alJ 
l2006i) . In the next section, we review the main properties 
of the general solution of eq. J28> for concentration domains 
of arbitrary geometry, which we subsequently specialize to 
the geodetic context by imposing the circular symmetry of 
the double-cap polar gap. 



3 SLEPIAN'S SPHERICAL PROBLEM 

In Slepian's problem, the concentration of a bandlimited 
function g given by 



9 ^^9lmYlm, 



9lr 



gYim dn, 



(29) 



to a region R of area A on the unit sphere fl is expressed as 
the norm ratio, eq. ( I28> . Maximization of this concentration 
criterion can be achieved in the spectral domain by solving 
the algebraic eigenvalue problem 



Dg = Ag, 



(30) 



where g is the {L + l)^-dimensional spherical harmonic co- 
efficient column vector 

g = (.900 • • • girn ■ ■ ■ 9LlV (31) 

and D is the (L + 1)^ x (L + 1) ^-dimensional matrix 

/ -DoO.OO • • • Dqq 



D 



(32) 



whose elements Dim I'm', < I < L and — / < m < /, are 



YlmYl'm' dn. 



(33) 



The obvious symmetry = D guarantees that the eigen- 
vectors g^, g2, . . . , g(i_|_i)2 are mutually orthogonal. We 
choose them to be orthonormal: 

sldp = 6ai3 and g^Dg^ = XaSap- (34) 

The resulting Slepian functions g are orthonormal over the 
whole sphere il and orthogonal over the region R: 

/ gagpdQ. = 5afi and / gag/s dfl ^ XaSa/s- (35) 
Jn Jr 
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The leftmost equations in eqs ( I34> -( I35> correspond to the 
conditions of eq. ( I23> and guarantee that the solution indeed 
forms a valid orthonormal basis. The rightmost equations il- 
lustrate the so-called double orthogonality of the Slepian ba- 
sis (Gilbert & Slepian 1977), which, as we will see in a later 
section, is a central feature of their utility for the geodetic 
estimation problem. 

An approach equivalent to the maximization of eq. j28t 
is to find broadband functions /i(f ) that are spacelimited to 
the domain R, but spectrally concentrated in a bandwidth 
interval Q <l < L. The concentration measure in this case. 



maximum, 



(36) 



2 



is satisfied by the eigenfunctions of a Fredholm integral 
eigenvalue equation in the spatial domain; 



R. 



(37) 



[ D{r,r')h{r')dn' ^ Xh{f), fe 
Jr 

The symmetric kernel of eq. ( I37l i depends only on the 
geodesic angular distance. A, between f and f ': 



^(f'f') = E(^~)^'(-f'). 



(38) 



The problems of finding bandlimited functions g concen- 
trated to a spatial interval or spacelimited functions h con- 
centrated in a spectral interval are completely equivalent. 
The domain of eq. ( I37> can be extended to the entire sphere 
in which case it applies to the bandlimited functions g: 



D{r,r')g{r')dn' ^ Xg{r), f e SI, 

Jr 



(39) 



We normalize such that the eigenfunctions g that maximize 
the spatial energy ratio ( I28t are identical, within the re- 
gion R, to the eigenfunctions h maximizing the spectral ra- 
tio (|36|i: 



h(r)^( .9(f) iff Gi? 
^ ' 10 otherwise. 

The relation 

L 

^Im ^ ^ I-^lm^l' m' QV m' 
I'm' 



(40) 



(41) 



expresses the coefficients hi„i, where < Z < cx), in terms 
of the coefficients gi„i, with < I < L. This is a straight- 
forward consequence of the definitions in eqs ( I29l l, i33i 
and ( I40> . and, by eq. i3Q\ . it amounts to him = ^9im when 
< I < L. The eigenvalues of eqs j30> or j37L 



1 > Ai > As • •• > A(i+i)2 > 0, 



(42) 



measure the quality of the spatiospectral concentration: the 
bandlimited function that is most concentrated inside R is 
gi, with Ai being the largest associated eigenvalue, and so 



on. The sum of the eigenvalues, or Shannon number, equal to 
the trace of D, defines a diagnostic area-bandwidth product 



iL+l] 
a=l 



[ D{r,r)dn = (L + l) 
Jr 



47r 



(43) 



Spherical Slepian functions of equal Shannon number are 
scaled versions of each other in the asymptotic limi t ^ ^ 
and L ^ oo with K held fixed dSimons et al.l200^ . When- 
ever the area A of the region i? is a small fraction of the area 
of the sphere, A <C 47r, i.e., when K {L + 1)^, there will 
be many more well excluded eigenfunctions with insignifi- 
cant eigenvalues (A « 0) than well concentrated eigenfunc- 
tions with significant eigenvalues (A « 1). If on the other 
hand, R covers most of the sphere so that A « 47r and 
K {L + 1)^, there will be many more well concentrated 
eigenfunctions than well excluded ones. 

The sum of the squares of the (L + 1)^ bandlimited eigen- 
functions g{f) is independent of position on the sphere: 



E Slir) = ^ 



1)^ 



a=l 



K 
~A' 



(44) 



Since the first K eigenfunctions gi,g2, ■ ■ ■ , gx have eigen- 
values near unity and lie mostly within R, and the remainder 
9K+1, gK+2, ■ ■ ■ , g(L+iy have eigenvalues near zero and 
lie mostly in the complementary region R ^ — R, the 
eigenvalue-weighted sum of squares is well approximated 
by 



{L+iy 

E 



K/A ifreR 
otherwise. 



(45) 



The terms with K < a < (L + 1)^ should be negligible, 
so it is immaterial whether they are included in the sum ( I45> 
or not. Taken together, the first K orthogonal eigenfunctions 
ga, a = 1,2, . . . , K, with significant eigenvalues Aq, « 1, 
provide an essentially uniform coverage of the region R. 
Rather than requiring (L + 1)^ basis functions to represent 
an arbitrary spatially concentrated bandlimited function, the 
first K = {L + l)^A/(47r) members of the Slepian basis 
provide a very reasonable approximation. 

We shall denote the operator localizing to the complemen- 
tary region R by D, its eigenfunctions by g and its eigenval- 
ues by A. It follows from the orthogonality relation (|6j that 
the elements of D are 



Im.l'' 



Su'd' 



(46a) 
(46b) 



The eigenfunctions of D are identical to those of D, but 
their ordering indices are reversed. The bandlimited func- 
tion that is most concentrated within R is most excluded 
from R, i.e. gi = .g(L+i)2, with an associated eigenvalue 



Ai = 1 - A 



, and so on. 



The localization operator D has an inverse satisfying 



E 

l"rn" 



D 



lm,l" 



l"rn" .I'tti' 



(47) 
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and for which D^^g = A^^g. For future reference, the in- 
verse of a weighted sum of the localization matrix and its 
complement obeys 



(D + 77D)-ig = [A + 7?(l-A)]- 



(48a) 



[Ac+rKl-Aa)] '5„/3, (48b) 



gI(D + 7yD)-ig,3 

for any weighting parameter -q. Finally, we may extend the 
validity of the definition (I33> to include all degrees I < oo 
and use eqs (I8t-(I10> to prove that 



Im.l'' 



' m" .I'm' 



Im.l' 



(49) 



4 AXISYMMETRIC DOMAINS 

In the previous section we showed that the optimally concen- 
trated bandlimited basis functions that are the solutions to 
the Slepian problem are found by numerical diagonalization 
of the operator in eq. i33\ . That this is in general possibl e 
for arbitrary geometries was shown bv lSimons et al.N2006h . 
However, the particular geometry of data acquisition on the 
sphere in the geodetic estimation problem (Figure^ allows 
for substantial simplifications of this general result. We dis- 
cuss the special case of finding concentrated basis functions 
on the latitudinal belt, the domain over which satellite mea- 
surements are made, via the concentration within the single 
and the double polar cap. As we have seen, the eigenfunc- 
tions on a domain R are identical to those on a complemen- 
tary spherical domain R, but with their ordering indices re- 
versed. Identifying R with the polar caps, rather than their 
complement, the belt, as we do - in this section and the one 
that follows only - greatly simphfies the equations. 



4.1 Concentration within an axisymmetric polar cap 

When the region of concentration is a circularly symmetric 
cap of colatitudinal radius 6, centered on the north pole, i.e. 



i? = {6* : < 6* < e}, (50) 

of area A = 27r(l — cos Q), the matrix elements of eq. ( I33> 
reduce to 



Im.l'' 



2-k5, 



XlrnXl'rn Siu dO . 



ram' I ^^Im^l'm 




(51) 



The Kronecker delta (5,„,„ renders the matrix D of eq. ( I32t 
block-diagonal. 



(52) 



D = diag (D^D^D^...,D^D^), 

where every submatrix D™ 7^ occurs twice due to the 
doublet degeneracy of ±m. Rather than solving the com- 
plete (L + 1)^ X (L + 1)^ eigenvalue equation ( 13 0> , we may 
solve a series of (L — m + 1) x (L — m + l) spectral-domain 
eigenvalue problems, one for each non-negative order m. 



where we have dropped the superscript identifying the order 
The eigenvalues belonging to every nonzero order, m > 0, 
occur twice. In eq. j53> the column vector g collects the 
spherical harmonic coefficients of order m: 

9 = {9rn ■■■ 91 ■■■ 9lY , (54) 

and the fixed-order matrix D is of the form 

D mm ' ' ' I^mL 

D = I : : I , (55) 

where, for a particular order < m < L, 

Ai'-27r/ Xi^Xv^sinede. (56) 

Various methods exist to ev a luate the elements o f eq. J56> 
jWieczorek & Simonsll2005t ISimons et alJl2006h . The im- 
portant point is that, while symmetric, and banded, the ma- 
trix D is never sparse. Its construction thus requires on the 
order of (L — m + 1)^/2 integrals each. 

We rank the L — m + l eigenvalues Ai, A2, . . . , Al-jh+i 
obtained by solving the fixed-order problem ( I53> so that 



1 > Ai > A2 > 



> A 



L — 771 + 



1 > 0, 



(57) 



and orthonormalize the eigenvectors gj, g2, • ■ • , gi-m+i 
as in eq. (I34> . The associated bandlimited eigenfunctions 

9i{^),92{d), ■ ■ ■,9L-m+i{d), are given by 



3 = E 9lXlm, 
l—m 



gi ==2tt / gXi,n sin 9 dO, 



and satisfy the colatitudinal orthogonality relations 

27r / ga9i3SinOde = Sap, 
Jo 

Jo 



(58) 

(59a) 
(59b) 



Dg = Ag, 



(53) 



The optimally concentrated spatial eigenfunctions g{r) for 
a given order —L < m < L are expressed in terms of the 
fixed-order colatitudinal eigenfunctions j58> by 

\/2 g{9) cosnKp if —L < m < 
9(0, (f>) = { g{6) if m = (60) 

\/2g{6)&\mrL(j) if 0<m<L. 

The fixed-order Shannon number 

L-m+l 

K,n^ J2 (61) 
a=l 

again is simply the trace of the fixed-order matrix D. 

We further note that the complementary fixed-order matri- 
ces are given by 

Du,=Su,-Du,. (62) 
The eigenfunctions of the fixed-order matrix D are identical 
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to those of D but appear in reverse order, and their eigenval- 
ues sum to one. The axisymmetric inversion formula analo- 
gous to eq. ( I47t is 



Y,D;i}Di„i,^5u', (63) 

/"=0 

and the axisymmetric analog to eq. (I49> is 

oo 



(64) 



I" —m 



4.2 Concentration within a double polar cap 

When the region of concentration is a pair of axisymmetric 
antipodal caps of colatitudinal radius Q, i.e., when 

R= {e -.Q <d <Q]\J {e -.-K -Q <e <-k], (65) 

of area A = 47r(l — cos 9), the reflection symmetry 

^im(^-e) = (66) 

checkers the fixed-order matrices D with zeroes, following 



Dn 



2tt 



l + (-l) 



(67) 



Comparison of eqs ( I56> and (|^} reveals that the eigenfunc- 
tions of the double-cap problem can be trivially obtained 
from the kernels belonging to the single polar cap. 

The spherical Slepian functions resulting from the diago- 
nalization of the double-cap kernel in eq. ( I67> are either even 
or odd across the equator Indexing their parity by p, as even 
{p = e) or odd {p = o), we modify eq. ( I58> to explicitly skip 
every other degree by using a primed summation symbol. 



9p 



giXim, 



where the lower limit irip is given by 
me = m and rrio = m + 1 , 
and the upper limit Lp is 

if m and L have the same parity 



L 

L - 1 



in 



opposite 

L — 1 if m and L have the same parity 
L opposite 



(68) 

(69) 

(70a) 
(70b) 



In this formalism, the coefficients that are required for are 
9m, gm+2, ■ ■ ■ , 9l if m and L are both even or both odd, 
and gm+2, ■ ■ ■ , 9l~i if itt- and L are of opposite parity. 
Likewise, the coefficients of go are gm+i, 9m+7i, ■ ■ ■ , 9l-i 
if m and L are both even or both odd, and gm+i, gm+3, ■ ■ ■ , 
gL if m and L have opposite parity. Eqs j66> and j68t then 
confirm that 

g,(d) = g,(7r - d) and = - 0). (71) 

While an equation of the form ( I53l l returns an alternation of 
even and odd functions with decreasing eigenvalues A, the 



indices of the matrix D may be permuted to form a block- 
diagonal form 

D = diag (De, Do) , (72) 
for which the half-size separate eigenvalue equations 

De9e = ^e9e ^nd DoQo = >^oSo (73) 

return exclusively even or odd solutions. This avoids round- 
off problems and speeds up the diagonalization. 

The slight perversity of our notation is that, in an all-even 
or all-odd approach as in eq. il3i . writing the fixed-order 
coefficient gi or indeed any expression involving the spheri- 
cal harmonic degree I, always has to be accompanied by the 



set of allowable degrees I, since I ~ m.p, 
depending on the parity. 



2,. 



5 THE MAGIC OF COMMUTATION 

While conceptually simple, the formalism presented in the 
previous section suffers from two important difficulties. 
First, assembling the matrices of eqs ( I56> and (|^} re- 
quires the calculation of Oil?) ma trix elements, by nu- 
merical integration or oth er means dWieczorek & Simons! 
I2OOI Isiinons et al.ll2()06l) . Second, and more importantly, 
when a large number of near-zero eigenvalues is present, e.g. 
when 0^0 and the complementary solutions are sought 
on R, the diagonaUza tion is rarely stable, as discussed by 
lAlbertella et al.l ( Il999h . In principle, any orthogonal set of 
solutions might suffice to solve the problem at hand, but 
those solutions will vary depending on the method of com- 
putation. The method outlined below always produces sta- 
ble, unique, solutions, and it does so at a speed which re- 
quires only 0{L) algebraic evaluations to construct the ker- 
nels. 



5.1 A commuting operator for the single polar cap 

In the case of the single symmetric polar cap, eq. (I37> can be 
rewritten as a series of fixed-order integral equations 

/ D{e,e')h{9')sm9' d9' = Xh{e), 0<9<e. (74) 
Jo 

each with an m-dependent, separable, symmetric kernel 

L 

D{e, 9') - 2^ ^ Xira{9)Xi„,{9'). (75) 



Building on the results derived bv lGilbert & Slepiaiil ([l^^, 
[Griinbaum et al. (1982) found a second-order differential 
operator that commutes with the convolutional integral op- 
erator of eq. (17 4> . For any < m < i, it is of the form 



T = (cos e - COS 9)Vl, + sin 9— - L{L + 2) cos 9, (76) 

d9 

where = d^/d9^ + cot9{d/d9) - m^{sm9)-^ is 
the fixed-order Laplace-Beltrami ope rator. The proof of th e 
commutation relation is sketched in ISimons et al 
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Since commuting operators have identical eigenfunctions, 
the spacelimited, fixed-order eigenfunctions h{9) can be 
found by solving the differential eigenvalue equation 

Th{0)=xh{e), 0<9<e, (77) 

where x 7^ A is the associated Griinbaum eigenvalue. 

Gr iinbaum's operator is a Sturm-Liouville operator 
dSim ons et al. 2006). Thus, eq. ( I77> has a simple and easily 
sorted spectrum, with an infinite number of distinct eigen- 
values Xi < X2 < ■ • ■ having an accumulation point at in- 
finity. The rank orderings of the eigenvalues xi,X2, ■ ■ ■ and 
the spatiospectral concentration factors Ai, A2, . . . , XL~m+i 
are reversed, so that the eigenfunction hi{6) associated with 
the numerically smallest eigenvalue xi, which has no nodes 
in the polar cap < 6* < 0, is the best concentrated fixed- 
order eigenfunction; h2{0), which has exactly one node, is 
the next best concentrated, and so on. 

Extending the domain of eq. Mil to the entire domain 
< 9 < n transforms the unknown functions from the 
spacelimited functions h again into the bandlimited func- 
tions g. Eq. Ml\ is then equivalent to the algebraic eigen- 
value equation 



Tg = xg, 



(78) 



where T is the (L — m + 1) x (L — m + 1) matrix with 
coefficients 



Tu 



2tt 



Xi^{TXi>,n)smdde. 



(79) 



Eqs (I78ll-(l79ll are completely equivalent to eqs ( I53> and ( I56> . 
Both matrices D and T are symmetric, D = D and T = T . 
In addition, they commute, DT = TD, so they have identical 
eigenvectors. In index notation. 



nl' 1 



(80) 



which can be used as a numerical check. 



There are a number of ways to evaluate the elements of the 
Griinbaum matr ix in eq . ( I79> , but the important result is that 
T is tridiagonal JSimons et al.i.20061) : 



Tu = -1(1 + l)cosQ, 
Tii+i = [lil + 2)-LiL + 2)] 



(2; + l)(2/ + 3)' 



Til' = otherwise. 



(81a) 
(81b) 



(81c) 



Eq. ( I78> can be used to find the {L — m + l)-dimensional 
eigenvectors g and thus the optimally concentrated polar cap 
eigenfunctions g{6) by numerical diagonalization of a tridi- 
agonal matrix T with analytically prescribed elements and 
a spectrum of eigenvalues x that is guaranteed to be regu- 
lar. Unlike the diagonalization of the original matrix D in 
eQ.( l53> . this procedure enables the stable computation of 
bandlimited functions that are optimally concentrated in a 
large rather than a small region of the unit sphere, as may be 
the case in geodesy. 



5.2 A commuting operator for the double polar cap 

Knowing that the solutions to the concentration problem for 
the double polar cap are either even or odd across the equa- 
tor, we may write the integral equation ( I37> . by analogy with 
eq. MA\ . as follows. Indicating the parity of the solutions by 
the subscript p, which takes the values p = e for the even 
solutions and p = o for the odd solutions, it can be seen that 



Dp{e,0') hp{e') sine' de' = xhp{0), 



(82) 



which is valid inside the double antipodal polar cap 

{e ■.o<e <e}u {e :tt -Q <e <Tr}, (83) 

and where the m-dependent kernel, analogous to eq. i75\ . is 



(84) 



As in eq. ( I68> . the primed summation skips every second 
entry, and the lower and upper limits are as in eqs dSgt-dTOt. 

Aga in basing ourselves on the results of lGilbert & Slepiaiil 
( 119771) and Griinbaum et al. ( 1982), we show in AppendixIXl 
that a Sturm-Liouville second-order differential operator 
that commutes with the convolutional integral operators of 
eq. (|82} is of the form 



(cos^ e - 

— Lp[Lp 



cos' 

+ 3) cos^ 9. 



2 cos 9 sin 9 



d9 



(85) 



The individual matrix operators Tp are once again tridiag- 
onal and symmetric and commute with the even or odd Dp 
of eq. (I73> . The bandwidth Lp is the same as in The 
elements of the double-cap Griinbaum matrices are 



i)cos^ e 
-2)(; + i) 

2 Sto^ 
~ 3 



2/ + 3 
- Lp(Lp 



-3)] 



T?, 



{2l + 3){2l - 1) 
Lp{Lp + 3)] 



1+2 



2/ + 3 




otherwise 



(86a) 



(86b) 



(86c) 



We again emphasize that, since we focus our attention sepa- 
rately on the kernels returning even or odd eigenfunctions 
or go, the degrees involved are restricted to 1,1' = nip, nip + 
2, . . . , Lp. Since every other degree in the matrix described 
by eq. ( I86> is skipped, both Tg and To are tridiagonal as in 
the single-cap case. As in eq. ( I73> we compute the even and 
odd eigenfunctions by separately solving 

Tpa„ = YrQ^ and T„( 



(87) 



Subsequently, we establish a single rank order of decreasing 
spatiospectral concentration: either per order, as in eq. (I57> . 
or across all orders, as in eq. ( I42> . 
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Figure 2. Colatitudinal dependence of the first six fixed-order, m = — > 4, eigenfunctions ga{6). o = 1 — » 6, bandfimited to L = 18, that ai'e 
well concentrated in the latitudinal belt extending ±60° on either side of the equator. The quality of the spatial concentration is expressed by the labeled 
eigenvalues Xa - None of the plotted functions show appreciable energy inside the complementai'y pair of antipodal polai' caps of radius © = 30° . 
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Figure 3. Colatitudinal dependence of the last six fixed-order, m = ^ 4, eigenfunctions ga {6), a = L — m + l^L — m — 4, bandlimited to L = 18. 
These are generally poorly concentrated in the latitudinal belt ±60° about the equator, except where the rank a exceeds the fixed-order Shannon number Km 
(examples in lower right). The functions that have the least energy inside of the equatorial belt, as shown by their low eigenvalues Xa, are best concentrated 
inside the complementary polar caps of colatitudinal radius © = 30° . 
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Figure 4. Three-dimensional spatial dependence of the first four fixed-order, m = — » 2, eigenfunctions ga{6), a = 1 — > 4, bandlimited to L = 18, well 
concentrated in the latitudinal belt extending ±60° on either side of the equator, as expressed by their eigenvalues Ac . Plot arrangement is as in Figure 111 
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Figure 5. Three-dimensional spatial dependence of the last four fixed-order, m = — > 2, eigenfunctions ga {9), o = Z/ — m + 1— >i — m — 2, bandlimited 
to L = 18, poorly concentrated in the belt ±60° about the equator, as expressed by their eigenvalues . Plot arrangement is as in Figurelsl 
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6 A SLEPIAN BASIS ON THE BELT 

Concentration within a single polar cap was treated ex- 
tensive ly by j Wieczpre k & Simons (2005) and Simons et alJ 
(I2OO6I) . We refer to their figures for illustrations. In this sec- 
tion we illustrate the solutions to the concentration problem 
when the concentration region contains all but an antipodal 
pair of polar caps of radius 0. Reverting to our notational 
convention in Section|2] we again use R to denote an equa- 
torial strip or latitudinal belt extending 7r/2 — north and 
south of the equator. Consequently, the antipodal pair of po- 
lar caps themselves is again defined to be the excluded re- 
gion R, in line with their role as the geodetic polar gap in 
which no satellite observations are available. 



6.1 Spatial-domain solutions 

The six eigenfunctions ga{6), a = 1 ^ 6 of fixed or- 
der < m < 4 that are most optimally concentrated in 
the latitudinal belt complementing a = 30° double po- 
lar cap are plotted in Figure |3 Their associated eigenval- 
ues Aq are listed to six-figure accuracy. The latitudinal belt 
ranges from 60° north to 60° south symmetrically about the 
equator. With the chosen bandwidth L = 18, the Shannon 
number defined in eq. i43i is K = {L + 1)^ cos Q « 313, 
which approximates the number of well concentrated eigen- 
functions with Awl. The best concentrated eigensolution 
of every order is a bell-shaped even function with no nodes 
in the belt. In keeping with the Sturm-Liouville character of 
the Griinbaum operator, every subsequent solution acquires 
one more node, so that the second best of every order is an 
odd function, the third is even, and so on. All of the eigen- 
values shown in Figure |2] calculated by numerically inte- 
grating eq. (I28> . are equal to one within six-figure accuracy, 
indicating that the concentration to the belt is nearly perfect, 
while both poles are almost completely excluded. Since the 
concentration region is very large, the calculation of these 
functions by any means other than the Griinbaum procedure 
described above will fail. 

With the parameters unchanged from Figure |2] Figure |3] 
shows the six worst concentrated eigenfunctions on the belt, 
ga{9), a = L — m + \ ^ L — m — A. These now naturally 
have almost all of their energy inside of the antipodal pair 
of polar caps of radius 9 = 30°. For the zonal functions 
of order m = 0, the even-odd alternation starting at q = 1 
with an even function in Figure |2]ends at a = L + 1 with an 
even function, since L itself is even. At m = 1, the sequence 
starts with an even function but ends at a = L with an odd 
function, at m = 2 with an even function at a = L — 1, and 
so on. Thus, the worst concentrated m = eigenfunction is 
even about the equator, the worst m = 1 eigenfunction is 
odd, and so on, in a pattern that alternates with increasing 
order. Had L itself been odd, the worst concentrated zonal 
function would have been odd, the worst m = 1 function 
even, and so on, reversing the pattern. 

Three-dimensional perspective views of the first four of 
the fixed-order m = — > 2 functions whose colatitudinal 
dependence we plotted in Figure|2are shown in Figure^] In 
accordance with eq. ( I60> the zonal m = Q eigenfunctions do 
not display any longitudinal zero crossings, since the num- 
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Figure 6. Cumulative energy of L = 18 bandlimited eigenfunctions 
concentrated inside of belts complementary to antipodal polar caps of 
radius = 5°, 10°, 15° and 20°. The Shannon numbers are K = 
360, 356, 349 and 339. The sums of squares gf{9, </>) + g^{9,(l>) + ■■■ 
and Ai(;^(6, (/)) + A2g|(6, (/))+ ■■■ are plotted versus colatitude 6, along a 
fixed arbitrary meridian 4>. Dashed lines show the full unweighted sums of 
(L + 1)^ terms, which attain the constant value K/A over the entire sphere 
0° < 6 < 180°. Solid lines show the eigenvalue-weighted partial sums of 
K terms and the full sums of (L + 1)^ terms, which are very nearly equal, 
and concentrated uniformly inside of the belt © < 9 < 180° — 0. 
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Figure 7. Cumulative energy of L = 18 bandlimited eigenfunctions con- 
centrated within circularly symmetric polar caps of colatitudinal radius 
= 5°, 10°, 15° and 20°. The Shannon numbers are if = 1, 5, 12, 22. 
The symbols used are identical to those of Figure l6l The solid lines show- 
ing the eigenvalue-weighted partial sums of K terms and the full sums of 
(L -I- 1)^ terms are very nearly equal, and concentrated uniformly within 
the pair of antipodal caps 0° < 6» < and 180°- & <9 < 180°. 



ber of longitudinal nodes follows the order m. Similarly, in 
Figure|5]we plot a three-dimensional rendering of twelve of 
the worst concentrated eigenfunctions of Figure|3] 

In Figure |6l we show the eigenvalue-weighted pointwise 
sums of squares J2a ^aQai^^^ 4>) for latitudinal belts comple- 
mentary to double polar caps of radii = 5°, 10°, 15°, 20°, 
with a bandwidth L = 18. The cumulative sums are con- 
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Figure 8. Eigenvalue spectra of the mixed-order L = 18 bandlimited operators concentrating within the equatorial belts that complement antipodal pairs of 
axisymmetric caps of radius O = 5° , 10° , 15° and 20° . The total number of eigenvalues is (L + 1)^ =361; only Ai — > Aio and A312 — » A361 are shown. 
Different symbols are used for the various orders — 18 < m < 18; juxtaposed identical symbols are ±m doublets. Top labels specify the rounded Shannon 
numbers K = 360, 356, 349 and 339. 



centrated inside of the latitudinal belt; solid lines in grey 
and black distinguish the sums carried up to the first K (the 
Shannon number) or all {L + 1)^ possible terms. In con- 
trast, the cumulative sums of the cap eigenfunctions, shown 
in Figure are concentrated within the double polar cap. 
The full unweighted sums ffaC^j 4>) of all (L + 1)^ terms 
(dashed black lines) are exactly K/A = [L + l)^/(47r) over 
the entire sphere in accordance with eq. ( I44t . and the expec- 
tation in eq. ( I45> is confirmed: inside of the concentration 
domain, the weighted sums approach K/A. 



der < TO < X are connected by lines, with each sequence 
offset horizontally by its order, and vertically by an arbitrary 
50 units, to faciUtate inspection. Thus, L + 1 eigenvalues 
Xi, X2, ■ • ■ , Xl+1 are plotted for to = 0, whereas a single 
eigenvalue xi is plotted for m = L. The spacing between 
adjacent fixed-order eigenvalues is roughly equant, without 
the numerically troublesome plateaus of nearly equal values 
apparent in Figure |8l This regularity is guaranteed by the 
Sturm-Liouville character of the Griinbaum operator Tp in 
eq. (|83. 



6.2 Eigenvalue spectra 

In Figure |8] we show the reordered, mixed-order eigen- 
value spectra for the concentration problem within the lat- 
itudinal belt between polar caps of colatitudinal radii 6 = 
5°, 10°, 15°, 20°. Once again the maximal spherical har- 
monic degree is L = 18. The rounded Shannon numbers 
K = 360, 356, 349, 339 lie in the middle of the steep, transi- 
tional part of the spectra, roughly separating the reasonably 
well concentrated eigensolutions (A > 0.5) from the more 
poorly concentrated ones (A < 0.5) in all four cases. There 
are many more functions that are well concentrated in the 
equatorial strip than there are that are concentrated inside of 
the double polar cap, as shown by the break at a = 10 in the 
abscissas. 

The corresponding Griinbaum eigenvalue spectra are 
shown in Figure |9] The ranked eigenvalues % for every or- 



6.3 Analytic continuation 

The Slepian functions 51(f), ... , g(i+i)2 (f) are defined on 
the surface of the unit sphere Q,. Together, they form a natural 
basis set for the expansion of potential fields and, in particu- 
lar, estimates of these fields, on a sphere of radius ||f || = 1, 
as in eqs. j20t-J21>. This new basis is localized: the support 
of the first Shannon number K basis functions lies mostly in 
the concentration region R, whereas the rernainder are con- 
centrated outside of this area of interest, in R. With satellite 
observations we are of course mostly interested in the sig- 
nal at some height above the surface of the unit sphere. We 
have previously derived an expression for the expansion of a 
field estimate at satellite altitude a, in eq. Mil . It is immedi- 
ately obvious from this equation that, even if we were only 
interested in the first K upward continued Slepian expansion 
coefficients of the estimate, we would still need to know and 
calculate the full set of {L + 1)^ coefficients at zero altitude. 
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Figure 9. Eigenvalue spectra of the fixed-order L = 18 bandlimited Griin- 
baum operators commuting witli tlie operators whose eigenvalues are shown 
in Figure Isl Separate sequences of eigenvalues Xii X2i • ■ • i XL—m+i for 
each angular order < m < L are connected by lines. Each sequence is 
offset horizontally by its order m, and vertically by 50 units per order. 



The full impact of this statement will not become clear until 
later in this paper, but to anticipate it we derive here a set 
of Slepian basis functions that are designed specifically to 
represent signals at an altitude. We can do this by interpret- 
ing the Slepian functions we have just constructed as poten- 
tial functions themselves. In that case their upward harmonic 
continuation onto a sphere larger radius ||r|j = 1 + a, where 
a > 0, yields new functions for which 



At = y^ffT ImYlm, 
Im 



g\hn = (1 + a) 



dlri 



(88) 



Had we instead defined the Slepian functions on the larger 
sphere to begin with, their analogues downward continued 
onto the unit sphere would be obtained as 



L 
Irn 



9lhn = (1 + a)'^^5z« 



(89) 



It is thus useful to define a symmetric, (L + 1)^ x (L + 1)^, 
downward continuation matrix A, whose elements are 



A 



Im.l'' 



= (l + a)'+i&.<5, 



(90) 



which allows us to restate the equations relating the upward 
and downward continued coefficients to each other concisely 
as: 

g|=A"^g, g = Ag|, (91a) 

91= Ag, g = A"V. (91b) 

The orthogonaUty relations of eqs ( I23t and ( I34t can be 
rewritten in terms of and in the form 



g|aA^g 



T/3 



gl ADAg.^ = XJ^p, (92a) 



gjcrt '^113= Sap, gXaA ^DA ^ Q'^is = (92b) 

and also 



and al 



(93) 



In this matrix notation we repeat eq. Ml\ as 

ST« = E (glA-'g/s) (94) 

[3=1 

We note for future reference that, although the transforma- 
tion matrix g^A^^g^j may be banded, it is not in general pos- 
sible to truncate it to circumvent having to calculate the full 
set of sp even if we are only interested in a truncated set of 
coefficients S|q. 

Finally, eq. ( I23t and eqs (|88}-(|93 can now be combined 
to prove the equivalent results: 



a— 1 ct— 1 a— 1 

which we use extensively in subsequent sections. 



(95) 



7 POTENTIAL FIELD ESTIMATION 

We return to solving the geodetic problem stated in Sec- 
tion |2] We are given noisy data, d, taken by a satellite at 
an altitude, a, over an incomplete sampling domain, R, and 
attempting to estimate the potential field, s, that gives rise to 
these observations, at its source level on the unit sphere, Vl. 
Although the source field has an infinite bandwidth, we will 
practically only be able to make bandlimited estimates of it, 
which we denote by J. The spectral limitation to the band- 
width L as well as the spatial restriction of the observation 
domain to the region R motivates our seeking an estimate 
in terms of a set of basis functions that are spatiospectrally 
concentrated, rather than using the non-localized spherical 
harmonics Yim of more conventional approaches. This new 
function set is the Slepian basis, ga, constructed in Sec- 
tions |3}|6l in a variety of geometries, but most notably for 
the axisymmetric case of a latitudinal belt around the equa- 
tor, and its complement the double polar cap, representative 
of the polar gap in geodesy. 

That the geodetic estimation problem is essentially a prob- 
lem of spatiospectral localization can be understood by 
considering a naive - and in practice unsuitable - estima- 
tion scheme. Suppose we construct estimate in the form of 
eq. (ED, 



(96) 



by minimizing its aggregate squared misfit with the data over 
the sphere, given by eq. M5\ . This amounts to solving the 
variational problem 



mmimum, 



(97) 



where the integration domain is the region R in which ob- 
servations are available. Substituting eqs ( I15l l and ( I25> into 
eq. (|^} and requiring the partial derivatives d^/dsim to 
vanish yields the condition 



s^YimdQ.^ I dYhndil, 

R JR. 



(98) 
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while the resuh d^^/dsf^-^ > as long as i? 7^ guarantees 
the convexity of the penalty function $. Inserting the repre- 
sentation into eq. (|98j and using the definition of 
the locaUzation kernel ( I33t and its inverse ( I47> . the estimate 
of the field coefficients at source level is given by: 



hm = (1 + a) 



I' rn' 



D 



imA'm' 



dYim' dO.. 



(99) 



Thus, the estimate depends on the inverse of the localization 
kernel D. It is therefore directly influenced by the size and 
the shape of the region of missing data, as well as by the cho- 
sen bandwidth. Since D tends to have a very low condition 
number (see, e.g., Figure|8}, finding a stable inverse D^^ is 
problematic: the geodetic inverse problem is ill-conditioned, 
as is widely advertised even without reference to the local- 
ization nature of the problem ( Rutl992akbi) . 

In the following sections we will derive alternative solu- 
tions whose qualit y we will judge using standard statisti- 
cal measures (e.g. ICox & HinklevllTgTl iBendat & Piersoll 
l2000h . The first will be the average of the squared differ- 
ence between a single estimate and the mean of all estimates 
over a set of realizations of the data, the estimation variance: 



v^{{s-{s)f) = {s^)-{sf 



(100) 



The angular brackets denote averaging over the ensemble 
of repeated observations, each observation being influenced 
by a different realization of the random noise. Similarly, we 
compute the difference between the mean of the estimators 
and the unknown signal, the estimation bias: 



b = 



(101) 



We refer to the difference between an estimate and the un- 
known signal as the estimation error: 



e = s — s. 



(102) 



Finally, we compute the sum of the variance and the squared 
bias term, known as the mean-square error, or mse: 



(103) 



For the moment we regard the unknown source signal s as 
the unique "truth", i.e. we consider s to be non-stochastic, 
although the data derived from it are contaminated by sto- 
chastic noise, see eqs ( I15t -( ll7t . 



8.1 Damped least-squares approach 

To stabilize the solution we amend the variational problem 
of eq. (|^} by including a weighted model norm: 

/ (s| - df dU + rj I S| dU = minimum, (104) 
J R J R 

where ?/ > is a damping parameter Retaining the spher- 
ical harmonic basis, once again we supply the bandlimited 
estimate 



Im 

and minimize ( I104t with respect to the unknown coefficients 
Sim- After minimal algebra, involving eqs ( I25t -( l26t . i33i 
and (|46|i-(|47|i, we obtain the spectral-domain solution, 

L 

Sim = (1 + a)'^^ ^ (Am, Cm' + ??Am,i'm') 



(105) 



X ; dYi>„i'dn, 

IR 



(106) 



which only holds at the degrees I < L, since, when / > L, 
no estimate is available, s/,„ = 0. The case where a = and 
?7 = 1, for which, from eq. i46i . D + D = I, the identity ma- 
trix, w as treated in some detail bv^Sneeuw & v an GeldererJ 
The integral over the data in eq. ( I106t is made ex- 
plicit by substituting eq. il5\ and using eqs ( I25t and ( I33t 
once again: 



dYim dCl 



'S]l' 



nYimdn. (107) 



Comparing eq. ( I106> to eq. ( I99l l, we now require the inverse 
of the weighted sum of the operator localizing to R and the 
complementary operator localizing to the region of missing 
data R. The addition of the small quantity 77D to the orig- 
inal matrix D improves its condition number. We postpone 
a discussion on determining the ideal value of the weight- 
ing parameter r/ but it is clear that the estimate of the field 
coefficients s/,„ in the form of eq. jl06t is now computable. 

In order to ascertain the statistical properties ( ll00t -( fT03t 
of the new estimate jl06t -( fTU7t we first calculate the aver- 
age of this estimate over all realizations of the noise. From 
eq. il6\ . this ensemble averaging of eqs ( ll06> -( fT07l i annihi- 
lates the random noise term, and we obtain 



8 SPHERICAL HARMONIC SOLUTION 

We have seen that a naive least-squares solution to the 
geodetic inverse problem in the spherical harmonic basis 
yields a solution ( I99> that is dependent on the inverse of 
the localization matrix and therefore in general impossible 
to stably compute. One of the many approaches to circum- 
vent this diffic ulty i s by adding a model norm to the penalty 
function (e.g. iHoe rl & Kennard 1970a b; Ma rauardlll 197(1 
lJacksonll979^ : eq. i91l only minimized the norm of the data 
misfit. In this section we discuss the solution to this so-called 
damped least-squares approach. 



(s;„ 



= (1 + a)'+l {Dim,l',n' + VDlrnMrn'Y 



I'm' 



X ^ Diim',l"m"S^l" 



(108) 



I" m" 



Again, the coefficients s;,„ are defined only in the degree 
range I < L. We note that, were the source signal to be 
similarly bandlimited, the coefficients sim obtained by un- 
damped (77 = 0) estimation would be equal to the true source 
coefficients s;,„. This follows directly from substituting the 
leftmost term of eq. ( I13> into eq. ( IIO8I 1 and using eqs MAI 
and ( I47> . The addition of the damping term (77 > 0) biases 
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the estimate away from the tr uth, hence the name "bias ed es- 
timation" for this procedure (iHoerl & Kennardll970lJ) . It is 
the price we pay to be able to calculate the estimate at all. 

There are other benefits as well. These are most easily 
seen by computing a spatial-domain representation of the es- 
timate using the Slepian basis, as in eq. \2\\ . Making use of 
the equivalence ( I95> , we write for the (bandlimited) estimate 



all realizations of the noise is given by 



(L+lY 

Q = l 



(109) 



noting that the upward continued coefficients s-\a in the 
Slepian basis are calculated according to eq. \21\ . and the 
downward continued Slepian basis functions according to 
eq. (I89> . A Slepian basis expansion of the (broadband) ob- 
servations, combining eqs (I15> and (I25> . is given by 



E 



Ylra + ' 



(110) 



lm> L 



This equation allows us to find an alternative expression for 
the data integral ( I107> . for which we also use eqs ( I22t . i33\ 
and the double orthogonality of the Slepian functions i35\ . 
namely 





J nga dil^ 







oo 

E 

/'m'>L 



-^hji^l'rn' S]' I'm' 



(111) 



Inserting eqs (|27}, ( fT06l l and (fTTTl into eq. ( fTol and usmg 
the expressions ( I41> and ( I48> yields the spatial-domain esti- 
mate of the field as 



{L+ir 



(112a) 




X ^AqS|q + J ngadfl 



We have introduced the symbol A* (77) for notational con- 
venience. In the absence of damping. A* (0) = A^;^, i.e. the 
inverse of the concentration eigenvalue. Here, too, the neces- 
sity of damping is readily apparent: as the eigenvalues of the 
concentration operator, Aq, become vanishingly small, their 
inverse grows explosively, inflating the noise term and the 
term containing the signal at the unmodeled degrees I > L, 
and rendering the stable computation of the estimate J112t 
impossible. Adding the damping factor is a useful way to 
prevent this. 

The zero mean of the stochastic noise, eq. guaran- 
tees that the ensemble average of the spatial estimate over 



{L+iy 

E 

Q = l 



(113) 



An 



'To + E ' 

lm>L 



'Till 



Eq. ( I113t can be combined with eqs ( I95l l and ( I20t to show 
that for bandlimited source fields and in the absence of 
damping, the mean of the spatial estimate (s(f )) is identical 
to the source field s(f ), even if the estimate s(f ) is impossi- 
ble to compute stably without the damping term. 

The introduction of the damping term stabilizes the solu- 
tion at the cost of added bias. Following eq. dlOU the latter 
is calculated by subtracting the full representation of the sig- 
nal i2Q\ from eq. (I113> . making use of eqs ( I48> and ( 195 > . The 
spatial estimation bias is then given by 



6(f) 



(L + lf 

(L+lf 00 

+ E ^a{v)9lair) ^ 

a—1 hn>L 



Im 



^ SlmYlTnir). 



(114) 



lm> L 



We have brought forward the damping parameter r/ by us- 
ing the identity A*A — 1 = — — A)A*. In the absence 
of damping [t] = 0), the first term in this equation vanishes, 
leaving us with the unavoidable broadband leakage (the sec- 
ond term) and bias due to making bandlimited estimates of 
broadband fields (the third term). 

An expression for the estimation variance from ( I100> is 
obtained by squaring eq. (II 12> and averaging the result, us- 
ing the properties of the noise ( I16> -( I17> an d eq. i35\ . and 
subtracting from the result the square of eq. ( I113> . The spa- 
tial estimation variance is 



viv) = N J2 UKM'9Ur). 



(115) 



We note that eq. dl 15> is the only one thus far to assume that 
the power spectrum of the noise is white, of magnitude N. 
And one more time the necessity of the damping is apparent: 
in its absence, the estimation variance strongly amplifies the 
measurement noise. At the price of introducing additional 
bias, damping prevents this. 



8.2 A bandlimited white stochastic source 

In the previous section we have derived expressions for the 
average esti mate o f the spherical harmonic field coefficients, 
(s;,„), in eq. ( I108> . and for the aver age o f spatial expansions 
of the estimated field, (s(f )), in eq. ( I113> . The averaging was 
over the different realizations of the stochastic noise process. 
Both expressions are valid in the most general sense; the 
only condition being that the average over all reahzations 
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of the noise, {n{r)), is zero. No further assumptions are nec- 
essary. We have drawn attention to the fact that without the 
damping term, both estimates are nearly impossible to cal- 
culate. However, in that case, they are unbiased when the 
source signal itself is strictly bandlimited to within a band- 
width L identical to that of the estimate. 

We can make this explicit by postulating that the geophys- 
ical signal expressed as eq. (I12t or eq. (I20> has spherical 
harmonic expansion coefficients that vanish outside of this 
bandwidth: 



for L <l < oo. 



(116) 



We will work with this contrived geophysical signal for the 
simple reason that no amount of sophistication can cure the 
fact that forming harmonically truncated estimates leads to 
multiple bias terms, as can be seen from eq. (II 14t . Under the 
condition (^TCJl, eq. jl08> becomes 

L 

{Shn) = (1 + a)'^^ X! {^hn,l'7n' + 'nDlm,l'ra') 
I'm' 

L 



X 



E 

l"m" 



Dl'rn' ,l"m"S]l"m" , 



(117) 



from which, using eqs (I14> and (I47> . we derive immedi- 
ately that the undamped estimate of the coefficients, given 
by eq. ( I99> . is unbiased: 



{kra) Shn if Tj = Q. 



(118) 



Similarly, using eq. ( I95> . eq. ( II 13> can be transformed under 
the same condition dl 16> into 



(119) 



(L+l)^ 



from which, with A* (0) = A^,^, the undamped spatial esti- 
mate of the field, given by eq. (I96> . is unbiased: 



(s(f )) s(f ) if = 0. 



(120) 



Indeed, under the condition dl 16> . the only term left in the 
bias equation (II 14> is directly, though not linearly, depen- 
dent on the damping term -q: it is 



(L + lf 

K^)^-V E (1 - ^")A;(?7)Sa5a(f). 



(121) 



Although we can calculate the mean-square estimation er- 
ror ( I103> exactly from eqs (II 15> and ( I121> . we will gain 
additional insight when we cease to consider the unknown 
signal as a non-stochastic signal. The source signal s(f ), un- 
til now, has been considered to be "given": we have simply 
assumed it is of the form (I12> and attempted to estimate its 
true unknown coefficients s;,„ from incomplete and noisy 
observations. All averaging in the construction of the bias 
and variance terms was carried out over the different realiza- 
tions of the noise n{r), which we took to be a white stochas- 
tic process. By now considering the geophysical signal, as 
well, to be a stochastic process, we shall calculate the mse 
after an additional round of averaging, this time over the var- 
ious realizations of s(f ), should they be available. Instead of 



eq. ( I103> we thus write 



(122) 



where the angular brackets now denote an average over the 
ensemble of signals. Strictly speaking we should write ((e^)) 
but we eschew the double brackets in the interest of nota- 
tional simplicity. 

We notice from eq. ( I121> that to compute (fe^) we shall re- 
quire the covariance (sqS^j) of the expansion coefficients of 
the field in the Slepian basis. To facilitate the treatment and 
for easy comparison with the assumed white power spec- 
trum of the noise process, we shall consider a bandlimited 
source signal that is "whitish", i.e. white within the band 
I < L, such that its covariances in the spherical harmonic 
and Slepian bases, respectively, are given by 



{simSl'm' ) 



S Su'Sjnm', 
S 



(123a) 
(123b) 



while noting that, as far as the spatial covariance of this sig- 
nal concerned, 

(s(f)s(f')) = SD{r,r') ^ SS{r,r'), (124) 

as can be deduced by combining eq. il2\ with eq. ( I123> and 
using eqs (|8|i-(|9j and ( I38> . A last assumption introduced 
here is that the noise is wholly uncorrelated with the signal: 

(s(f)ji(f')) = 0. (125) 

The average of the squared bias term ( I121> under these ide- 
alized assumptions is 



(L+iy 

E 



(1 



(126) 



and the mean-square estimation error, following eq. ( I122> . 
is formed by combining this result with the expression for 
the variance in eq. dl 15> . The latter expression is unchanged 
even if the source signal is stochastic, as long as the noise 
is uncorrelated with the signal, eq. dl25> . Thus, the mean- 
square error of the bandlimited estimation of a bandlimited 
white source field from incomplete observations at an alti- 
tude in the presence of white noise is given by 



(127) 



(62(f)) = N J2 K[Kiv)?9U^) 

a = l 

{L + lf 



All (L + 1)2 basis functions are required to form the mse. 
The first term in the expression for the mse is the variance: 
it is the only term that depends on the noise. We have seen 
that without damping (77 = 0) this term becomes unman- 
ageably large: the addition of damping counteracts this. In 
addition, the estimation variance also varies with the obser- 
vation height a above the unit sphere: as a grows, so do the 
downward continued Slepian basis functions (f ), and with 
them, the noise. The second term in the mse is due to bias. 
This is the only term that depends on the characteristics of 
the signal. It is independent of the satellite altitude at which 
the measurements are taken. 
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Figure 10. Spatially averaged mean-square en'or (as a percentage of the average mean-square signal strength) of the damped-least-squares spherical harmonic 
solution to the geodetic estimation problem for a bandlimited white signal and white noise. The Q-average mse (black curves and left ordinate) is the average 
over the entire sphere; the R-average mse (grey curves and right ordinate) is the average over the region of observation, the equatorial belt complementary 
to the polar gap of radius = 10°. The bandwidth of the signal and its estimate is L = 45. The measurement altitude is a = 0. The signal-to-noise levels 
shown are S/N = 4, 6, 8 and 10. We plot the normalized average mse values as a function of the damping parameter r], and indicate by vertical lines the 
values of rj that minimize them. The range of the i?-average is much reduced compared to the ^-average values. Both ordinates are truncated below at N/ S, 
the mse value when the observation region is the entire sphere, R = Q. 



8.3 Optimal damping level 

To illustrate the behavior of the mse in eq. (I127> we will 
focus on the case where the measurement altitude is a = 0, 
hence ga — Oia ~ 9^ a- This simplifies the expressions to; 

(L + lf 

(e2(f)) = J2 (128a) 
T^aiv) = [A;(?/)]2 [iVA,+ 7725(1 -A„)2] . (128b) 



The function TZairj) combines the effects of data noise, 
damping, signal strength, and measurement geometry. We 
will compare the mean-square error with the mean-square 
signal strength over all realizations, which is given by 



An 



(129) 



The result ( I129> is obtained by combining eq. ( I124t with 
the definition j38t at f f — 1. We calculate the following 
two quantities. First, a normalized spatial average of the mse 
given by the ratio of the mean square error ( I128> to the mean 
square signal strength ( I129> , both averaged over the entire 
sphere H. Using the orthogonality conditions ( I35> this "17- 
average mse" is given by 



S{L 



(130) 



Second, a scaled "i?-average mse" is given by the ratio of the 
same quantities, averaged over the covered region R. Using 
eq. ( I35t and the definition ( I43t of the Shannon number K, it 



IS 



(131) 



■{r))dn 



Both quantities are shown in Figure^! for a double-cap po- 
lar gap of 6 = 10° and a bandwidth L = 45. They are 
plotted in different panels for different signal-to-noise ratios 
S/N = 4, 6, 8 and 10 as functions of the damping parameter 
7] = ^ 1. We show eq. jl30> in black, with the scale on the 
left of the panels, and eq. ( I131> in grey, with the scale on the 
right hand side. The range of 0-average mse values shown is 
much larger (5% in all four panels) than the equivalent range 
in i?-average mse values (0.4% in all panels): the effects of 
damping on the overall mse over the entire globe are much 
more pronounced than its effects on the mse averaged over 
the region in which data were collected. The ordinate is trun- 
cated to aid the visualization. The maximum i?-average mse 
is {N/ S) I cos 6 which is attained when 77 = 0. This can be 
verified by noting that 7?.q(0) = , using the definition 

of the Shannon number ( I43> . and noting that the area of the 
covered region is equal to A = 47r cos 9. Thus, at a given 
signal-to-noise ratio only the size of the polar gap controls 
the upper bound on the i?-average mse. A lower bound for 
all damping levels is found at full coverage, i? = il. It thus 
applies to both measures of the average mse. Indeed with- 
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Figure 11. Mean-square error (as a percentage of the mean-square signal strength) of the damped-least-squares spherical harmonic solution to the geodetic 
estimation problem for a bandlimited white signal and white noise. As in Figure [Tol the signal-to-noise levels shown are S/N = 4, 6, 8 and 10, the mea- 
surement altitude a = 0, the polar gap consists of caps with radius © = 10°, and the bandwidth of the signal and its estimate is L = 45. As a function of 
colatitude, for an arbitrary longitude, we plot the mse of the undamped solution (r; = 0), of a heavily damped solution (r? = 1), and of the solution at the 
damping level which minimizes the normalized average mse over the unit sphere, i.e. for the values t^q = 0.37, 0.22, 0.15, and 0.12, that were marked by 
black vertical lines in Figure fTol The mse is symmetric about the equator. The ordinate is truncated at 100%; the mse of the undamped solution in the region 
of the polar gap exceeds this value by several orders of magnitude. 



out a polar gap, 6 = 0, K = (i + if, A;(7?) = A„ = 1, 
T^a iv) = the scaled average mse curves never drop 

below N/S, which we use as a lower cutoff for the vertical 
axes. 

A statistically desirable e stimator (e.g. ICox & HinklevI 
ll974HBendat & PiersoB200Clh is one that is unbiased and ef- 
ficient, i.e. it minimizes the mean square estimation error 
We have seen that sacrificing the unbiasedness by introduc- 
ing damping removes the obstacles in computing the esti- 
mate in the first place, and reduces the estimation variance. 
We can calculate the damping level that is overall optimal 
by minimizing the mse jl28> with respect to the damping 
parameter rj. However, minimization of the i?-average and 
f2-average mse will yield slightly different optima. Minimiz- 
ing, eq. jl30> , the normalized mse over the entire sphere we 
obtain an optimal damping coefficient Tyjj given by 



^ EiiV^'[A*(r?oFA„(l-A„)2- 



(132) 



Likewise, minimization of eq. ( I131> . the normalized mse 
over the region of coverage, yields an optimal damping co- 
efficient ?7fl given by 



N EiiV'lA;M^A^(l-A.) 
^ ELiV^'[A*(^fiFA2(l-A.)2- 



(133) 



Although the unknown optimal damping levels ijfi and rjj^ 
appear on both sides of eqs ( I132> and ( I133t . their values can 



be easily computed by iteration. They depend on the mea- 
surement geometry, the damping, and the signal-to-noise ra- 
tio. In Figure^l rjn and rjj^ are shown as black and grey ver- 
tical lines, respectively. At high signal-to-noise ratios both 
can be approximated as 77^7 = 77^ ~ N/S <C 1. 

When the coverage region is axisymmetric the mse jl28t 
is independent of the longitude, as can be deduced from 
eq. (|60|l. Thus, in FigureFTTIwe plot (e^(6')) /{s^iO)), in per- 
cent, for different signal-to-noise ratios, as a function of co- 
latitude and for various damping levels: i.e. in the undamped 
(77 = 0), fully damped (77 = 1) and optimally damped case 
ill = '7si)- The vertical axes are truncated at 100% as the 
undamped values exceed this value by many orders of mag- 
nitude. 



9 SLEPIAN BASIS SOLUTION 

In the previous section we expanded the estimate of the sig- 
nal into a bandlimited spherical harmonic basis and per- 
formed a damped least-squares inversion for the unknown 
coefficients. This estimation procedure resulted in a biased 
estimate, but the damping prevented the detrimental amplifi- 
cation of the measurement noise. We derived expressions for 
the optimal level of damping required for "whitish" signals 
measured at zero altitude. Adding a small amount of bias 
made the estimate computable and reduced its variance. We 
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used the (downward continued) Slepian basis to find expres- 
sions for the resultant damped spherical harmonic estimate 
in the spatial domain and to find its bias, variance, and mse. 
Using the Slepian basis greatly simplified the expressions 
because of the fact that, as opposed to the spherical harmon- 
ics, the Slepian functions are orthogonal over both the entire 
sphere and the closed subdomains over which, by design, 
their energy is optimally concentrated. 

Alternatively, we might have sought an estimate that is 
expressed in the spherical Slepian basis at the start. As we 
have seen, the first K Slepian eigenfunctions, where K is 
the Shannon number ( I43> . provide an excellent coverage of 
the region of observation. This implies that their associated 
eigenvalues A are close to unity, avoiding any problems with 
their inversion. In this section we will explore the effect on 
the geodetic solution of using a truncated Slepian basis, con- 
sisting of the J basis functions that are best concentrated 
over the region of satellite observation. Even if J = K ap- 
pears to be a natural choice, we will determine the trunca- 
tion level J by optimization of the mean-square estimation 
error, as we did to find the optimal damping parameter in the 
damped least-squares spherical harmonic approach. 



9.1 Truncated Slepian function approach 

The original undamped problem posed in eq. ( I97L 

/ (s| - df dfl = minimum, (134) 
Jr 

is now solved by expanding the estimate in the downward 
continued truncated Slepian basis 



9lc 



(135) 



and minimizing eq. with respect to the estimation coef- 
ficient S|cf The second derivative of eq. (I134> is always pos- 
itive. After minimal algebra, using eq. (125 > and the double 
orthogonality ( I35> . the expansion coefficients in eq. ( I135> 
are obtained from 



JR 



dn. 



(136) 



This result can alternatively be derived by substituting 
eq. ( I106> of the damped spherical harmonic approach into 
eq. (I27> . setting r/ = 0, and using eqs (I22> and ( I48> . We pur- 
posely chose an estimate the form ( I135> to find the truncated 
expansion coefficients S|q, in their upward continued form 
and multiplying the downward continued Slepian functions 
Qla, rather than simply expressing eq. ( I99> in the Slepian 
basis Qa- In the latter case, as can be readily verified by 
combining eq. ( I99> with eqs (I24> . (I22> and j48> . every one 
of the expansion coefficients Sa would depend on a linear 
combination of all {L + 1)^ terms A^^ dg/^ dfl through 

a matrix term g^Ag^ whose kind we have encountered in 
eq. ( 194k This would therefore invalidate the method of trun- 
cation as a means to avoid the difficult-to-compute and un- 
necessarily influential large inverse eigenvalues. By choos- 
ing the representation ( I135> instead, we take advantage of 
eq. ( I95> to juxtapose upward and downward continuation. 



(1 + a)'+^(l + a)^'^^, thereby canceling their effect alto- 
gether; eq. ( I136> shows that every coefficient S|q only de- 
pends on the inverse eigenvalue at the same rank a. The ef- 
fect of the measurement at altitude has not disappeared: it is 
now contained in eq. ( I135> in the basis gia, of which only 
the first J functions are required. These are calculated via 
eq. ( I89> and ultimately, by the stable Griinbaum algorithm 
central to our analysis. 

The data integral ( I136> can be calculated by substituting 
into it eqs ( fTTOl . (|35}, (EJJ, (|33} and (iD, to yield 

/ dgadn ^ XaS^a + / ngadVt+ ha.lmS^lm- 
JR JR r 



lm>L 



(137) 



Averaging the expressions ( ll36> -( fT37l over many estimates 
annihilates the influence of the random noise by virtue of 
eq. ( I16> . and gives 



lm>L 



(138) 



Combining eqs ([03}-<[^3 yields the estimate in the space 
domain, 

,/ 

s(f) = ^A-ig^„(f) (139) 

a=l 

X AaS|Q + / ngadn+ ^ ha.hnS^lm , 
V lm>L J 

which, reassuringly, amounts to the truncated but undamped 
(r/ = 0) version of eq. dl 12> . As before we can eliminate the 
noise term by averaging over many realizations, to obtain 

(s(f)) = ^.9i„(f) 

X + ^ ha^lmS-\lrn ■ (140) 

\ lm>L I 

The estimation bias, following eq. (I101> . is obtained by sub- 
tracting from eq. ( I140> the representation of the signal ( I20> 
and using the equivalence ( I95> : 

{L + \f oo 

b{r) = - ^ Saga{r) - ^ Simi^im(f) 
a>J lm>L 

J OQ 

+ ^K^9la{r) ^ ha^lmS^lm- (141) 
Q=l lm>L 

Without truncation of the Slepian basis function set, i.e. 
when J = (L + 1)^, the first term in this equation van- 
ishes. The remaining contributions arise due to forming 
bandlimited estimates of broadband signals, leading to un- 
avoidable leakage and broadband bias. Comparing eqs dl 14> 
and ( I141> we discover the parallel roles of damping and trun- 
cation. The introduction of the damping parameter rj adds 
an extra bias term to the expression dl 14> . and reduces the 
size of the leakage term by which the coefficients ha^im, 
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Im > L of eq. i36\ make the influence of the signal out- 
side the bandwidth felt, but it is powerless against the bias 
due to the bandlimited approximation of the broadband sig- 
nal, which is simply that portion of the signal that is outside 
the bandwidth L. Similarly, increasing the Slepian trunca- 
tion level by the reduction of J from (L + 1)^ in eq. 
introduces a new term in the expression for the estimation 
bias, and reduces the effect of the leakage term containing 
the coefficients h, but it is again no match for the remaining 
broadband bias from the bandlimitation of the estimate. 

An expression for the estimation var iance , eq. (llOOt . is ob- 
tained by squaring and averaging eq. ( I139> . using the noise 
properties ( I16> - (I17> and the orthogonality of the Slepian ba- 
sis functions i35\ . and subtracting the square of ( 114m . The 
resulting variance is 



(142) 



a=l 



This expression is again the first in this section in which we 
have used the white noise assumption, and once again it will 
be valid even if the source signal is considered stochastic as 
long as eq. ( I125> holds. Comparison of the variance expres- 
sion in this truncated Slepian basis approach with eq. il 15> 
obtained via the damped spherical harmonics method val- 
idates our approach. Without damping, when 77 = in 
eq. ( II 15> , or without truncation, J = (L + 1)^ in eq. ( I142> . 
both expressions are identical. Much like the damping term, 
the truncation of the basis set to its first J < {L + 1)^ 
elements reduces the estimation variance by checking the 
growth of the terms A^^. The more severe the truncation, 
the lower J, and the lower the variance becomes. 



9.2 A bandlimited white stochastic source 

Once again, we now focus on geophysical signals that are 
white within a bandwidth L as expressed by eqs (II 16> 
and ( I123> . This assumption transforms eq. ( I138I I into 

(S|a) = (143) 

illustrating the fact that an estimate of the form il 'ibt is spec- 
trally unbiased. Just as our analysis of the damped spherical 
harmonic method showed that for bandlimited source fields, 
the undamped, i.e. rj = 0, estimate of eq. ( I99> is incom- 
putable due to the ill-conditioning of D^^, but unbiased, as 
shown by eq. ( I118> . we have now shown that the untrun- 
cated, i.e. a = 1 ^ {L + 1)^, Slepian basis estimate of 
eq. ( I136t is incomputable due to the growth of the eigenval- 
ues A^^, although it, too, is unbiased, as shown by eq. ( I143> . 
The damping term makes the estimate computable but bi- 
ased, just as the truncation of the eigenvalues prevents the 
blow-up of their inverse at the cost of added bias. 

In the spatial domain, using eqs (^^J and ( I95> , the average 
over all estimates is then 



(144) 



In the absence of truncation, J = {L + 1)^, the spatial esti- 
mate of the form ( I135l l is similarly unbiased: 

(s(f)) = s{r), (145) 



which we may again compare to the unbiasedness ( I120> of 
the undamped estimate ( I96> . Explicitly, under the condition 
dl 16> , the only contributing term in eq. ( I141l i is given by 

(L+l) = 

Kf) = - sc.5a(f). (146) 

This term decreases with increasing J, and vanishes alto- 
gether when J = (L + 1)^. It can be compared to eq. ( I121> . 
The average over all realizations of the signal of the squared 
bias, for a "whitish" signal with covariance ( I123> , is given 
by 



{L+iy 

a>J 



{b\r)) = S 



.9a (f), 



(147) 



which should be compared with the corresponding eq. ( I126t 
in t he d amped spherical harmonic case. From this and 
eq. ( I142t we can calculate the mean-square estimation error 
following eq. ( I122t . which is now 



.7 

E 



(L + l) = 



(148) 



a>J 



The mse of the untruncated Slepian basis approach 
and that of the undamped spherical harmonic estimation 
method ( I127l i are identical. This of course is a direct conse- 
quence of the fact that both bases are related to each other by 
the orthonormal transformation eqs ( I22t -( l23t . Of note is the 
very different form of the damped and truncated expressions, 
eqs ( I127> and ( I148t . for the mse. Whereas eq. ( I127> consists 
of a weighted sum of all basis functions Q! = 1^(L + 1)^ 
in a manner that appears to mix the influence of the noise, 
the damping, and the signal, the truncated expression ( I148t 
has disentangled the effects of the noise and the signal by 
distributing the influence of the variance over the basis func- 
tions a = 1 ^ J that are well concentrated inside the 
measurement area, and the effect of the bias over those 
a — J+1-^{L + 1)^ that are concentrated in the re- 
gion of missing data. To the one piece that is missing, the 
decision on where to truncate the data by the choice of J, 
we now turn. 



9.3 Optimal truncation level 

To illustrate the behavior of the mse in eq. ( I148> we again 
focus on the zero-altitude case, for which 



(L+l) 



(149) 



a>J 



In order to find the optimal truncation level, we consider 
the full-sphere and coverage-domain average mse (normal- 
ized by the corresponding quadratic signal averages) as in 
the damped spherical harmonic approach. Using eqs ( I129> 
and ( I35> we find from eq. ( I149> that the f2-average mse in 
the truncated Slepian case is 
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Figure 12. Spatially averaged mean-squai'e en'or (as a percentage of the average mean-square signal strength), of the truncated Slepian function solution to 
the geodetic estimation problem for a bandlimited white signal and white noise. The f2-average mse (black curves and left ordinate) is the average over the 
entire sphere; the i?-average mse (grey curves and right ordinate) is the average over the region of observation, the equatorial belt complementary to the polar 
gap of radius = 10°. The bandwidth of the signal and its estimate is L = 45. The measurement altitude is a = 0. The signal-to-noise levels shown are 
iS/AT = 4, 6, 8 and 10. We plot the average mse values as a function of the truncation rank J, and indicate by sohd vertical lines the values Jq = J/j that 
minimize them, and the Shannon number K = 2084 by the dotted vertical line. The range of the i?-average is much reduced compared to the Q-average 
values. Both ordinates are truncated below at N/S, the value when the observation region is the entire sphere, R = Q. The abscissa shows truncation levels 
from J = ii" — 50 up to the untruncated case, J = (L -I- 1)2. In that case the mse values are identical to those of the undamped (rj = 0) spherical harmonic 
case shown in Figure lTol 



Using the definition ( I43> of the Shannon number, we Hke- 
wise find the i?-average mse: 



N J 1 

Ik ^ K 



(151) 



a>,J 



Both quantities are plotted in Figure^jfor different signal- 
to-noise ratios S/N = 4, 6, 8 and 10 and with the other pa- 
rameters unchanged from those of Figure^] a double-cap 
polar gap of O = 10°, a bandwidth L = 45. In black, with 
the scale on the left of the panel, we show eq. ( llSOI l as a 
function of the truncation level ranging over J = A' — 50 ^ 
(L -I- 1)^. The abscissa is inverted since J = {L + 1)^ cor- 
responds to a situation without Slepian truncation, and as 
J decreases, the degree of truncation increases. In grey, we 
plot eq. J151> , with a much reduced scale on the right hand 
side. The range of fi-average mse values shown is much 
larger (5% in all four panels) than the equivalent range in 
i?-average mse values (which varies from panel to panel but 
is less than 0.8%): the effects of truncation on the overall 
mse over the entire globe are much more outspoken than its 
effects on the mse averaged over the region in which data 
were collected. This behavior mimics the one seen in Fig- 
ure^Jfor the damped spherical harmonic case. The ordinate 
is again truncated for clarity. The value of the untruncated R- 
average mse, attained when J = (i -I- 1)^, is (TV/ S) / cos Q. 
This follows from the definition of the Shannon number i43i 



and the area of the covered region, A = 4tt cos Q, and is 
identical to the corresponding value in the undamped spher- 
ical harmonic case. A lower bound is found at full coverage, 
when R = n,e^O°,K^{L + if and A„ = 1. In that 
case the minimal scaled mse, attained when J = {L + 1)^, 
equals N/S, as it does in the damped spherical harmonic 
case. We use this value as a lower cutoff of the vertical axis 
on the left and the right. 

We may obtain the truncation level that minimizes the il- 
average mse by minimizing eq. ( I150> with respect to ,/. This 
will yield a n opt imal truncation level Jq. Likewise, mini- 
mizing eq. ( I151> returns the truncation value Jfi at which 
the i?-average mse is minimal. Both minimization problems 
result in identical constraints on the eigenvalue of the J-th 
eigenfunction beyond which we truncate: 



X, 



(152) 



which is implicit but solvable. In Figure the values 
Jn = Jr. identified by the top labels, are shown as a single 
solid black vertical line; the Shannon number, K = 2048, is 
shown by the dotted black line and the bottom labels. 

The mse jl49t for axisymmetric coverage regions is inde- 
pendent of the longitude, as can be understood from eq. ( I60> . 
Thus, in Figure^jwe plot {e^{0)) /{s^{9)), in percent, as a 
function of colatitude for various truncation levels: the un- 
truncated, ,/ = {L + 1)^, and optimally truncated cases, 
J = Jn = Jr, and the case truncated at the Shannon num- 
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Figure 13. Mean-square en'or (as a percentage of the mean-square signal strength) of the truncated Slepian function solution to the geodetic estimation problem 
for a white signal and white noise. As in Figure [Til the signal-to-noise levels shown are S/N = 4, 6, 8 and 10. the measurement altitude a = 0, the polar 
gap consists of caps with radius = 10°, and the bandwidth of the signal and its estimate is L = 45. As a function of colatitude, for an arbitrary longitude, 
we plot the mse of the untruncated solution, when J = {L -\- 1)^, of truncation at the Shannon number, J = K, and of the truncation levels which minimize 
the mse, i.e. = Jji = 2092, 2096, 2096, and 2098, that were marked by black vertical lines in Figure [ill The mse is symmetric about the equator The 
ordinate is truncated at 100%; the mse of the untruncated solution in the region of the polar gap exceeds this value by several orders of magnitude. 
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Figure 14. Mean square error (mse), variance and bias for the damped least-squares solution and the truncated Slepian approach. The antipodal polar gap has 
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ber J ~ K. The vertical axes are truncated at 100% since 
the untruncated values exceed this value by many orders of 
magnitude. 

Finally, in Figure O we plot the relative contributions of 
variance and bias for both the damped spherical harmonic 
and the truncated Slepian case. The top panels show the 
mean square error, the variance and the squared bias, ac- 
cording to the relation (I122> . as a function of the damping 
level (top left) or the truncation rank (top right). The bot- 
tom panels show the breakdown of mse, variance and bias 
in the spatial domain. For both estimation methods, the bias 
is predominantly concentrated in the areas over which no 
measurements are available, where it is generated by the 
power of the "missing signal". The variance, on the other 
hand, arises in the areas of coverage and is influenced by 
the power of the noise. Both estimation methods show very 
similar results. Over the covered area, the mse is nearly iden- 
tical, and in the uncovered region, the mse of the truncated 
Slepian case approaches that of the damped spherical har- 
monic case, but it is slightly higher It is remarkable that 
eqs (I128> and ( I149> . despite their different form, both give 
rise to a nearly complete spatial separation of bias and vari- 
ance. Only in eq. jl49> is this separation immediately obvi- 
ous by inspection: signal strength, and thus bias, affect the 
low-ranking Slepian functions, whose power is mostly con- 
centrated inside of the polar gap, whereas noise, and thus 
variance, affect the high-ranking Slepian functions whose 
power is localized to the area of satellite coverage. 



10 CONCLUSIONS 

Spherical Slepian functions provide a natural solution to 
the geodetic problem of having a polar gap in the satel- 
lite coverage of planetary gravitational or magnetic fields. 
Indeed, the ill-posed geodetic estimation problem of find- 
ing the source-level potential from noisy observations taken 
at an altitude over an incomplete region of coverage has 
natural connections to Slepian's spherical problem of spa- 
tiospectral localization. We have proposed a new method 
that expands the source field in terms of a truncated basis 
set of spherical Slepian functions, and compared its statis- 
tical performance with the traditional damped least-squares 
method in the spherical harmonic basis. The optimally trun- 
cated Slepian method performs nearly as well as the opti- 
mally damped spherical harmonic method, but it has the sig- 
nificant advantage of an intuitive separation of the estima- 
tion bias and variance over those Slepian functions sensitive 
to the uncovered and covered regions, respectively. The con- 
struction of Slepian functions over axisymmetric domains 
such as the latitudinal belt or its complement, the polar gap, 
previously dismissed as computationally unstable, has been 
shown to be eminently tractable. We have shown that the op- 
erator that bandlimits a field on the unit sphere and projects 
it onto the polar caps commutes with a Sturm-Liouville op- 
erator Its eigenfunctions, the Slepian functions, can be com- 
puted extremely accurately and efficiently by diagonalizing 
a tridiagonal matrix with analytically prescribed elements. 
The gains in ease, speed, and accuracy thus achieved makes 
the use of spherical Slepian functions in earth and planetary 
geodesy practical, as our examples have shown. 
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To confirm commutativity we are required to show that 



Dp{fj.,fi')Tphp{fi')dfj.' 



f TpDpifi, fi') hpifi') dn' . (A.6) 

Jb 



We first show that the left side of eq. ( IA.6t can be rewritten 

as 



APPENDIX A: GRUNBAUM COMMUTATION 

In this section we prove that the differential operator of 
eq. ( I85> . rewritten for i^i = cos 6 and b = cos Q, 

Tp = {b^~^l^)^l,-2^i{l-^?)^-Lp{Lp + i)^i'' 



d 

dfi 



dfi 



Lp{Lp + 3)fi^ 



(A.l) 



commutes with the integral operator acting on hp{fi') in 
eq. (|82}, 



(A.l) 



whose symmetric kernel, Dp{iJL, ^jl'), is given by eq. (I84> . 
rewritten here for convenience in the form 



Dp{fl,^i') = ClrnPlm{p-)PlmifJ-')- 



(A.3) 



The domain of eq. iA.2i is the interval ( l83t of the double cap 



{/X : 6 < ^ < 1} U {/i : -1 < ^ < -b}. 



(A.4) 



As in eqs (E}-(|5), Pim is the associated Legendre polyno- 
mial of degree I and order m, and Cim is a normalization 
constant. We remind the reader of our notation: hp{ji) is a 
colatitudinally dependent function that is limited in space to 
the antipodal polar caps of radius = cos^^ b. It is either 
odd or even about the equator, as indicated by the subscript 
p: 



hpin) = 0, -b> n>b, 



(A.5a) 
(A.5b) 
(A.5c) 



The solutions to eq. (IA.2> are functions hp{^) that are spec- 
trally concentrated in a spherical harmonic degree interval 
< / < L; the eigenvalue A is the quadratic measure of 
this concentration ( I36> . The primed summation symbol skips 
every other term in the interval from nip to Lp, which are 
both of the same parity, either even or odd. Depending on 
the requested order m and concentration bandwidth L of the 
solutions, TOp is either m or to + 1, and Lp is either L or 
L — 1, following eq. ( I70> . We further distinguish T acting 
on ji from T' that acts on ji' . 



Dp{^i,^i')r;hp{^l')d^Ji' 



I T;Dp{fi,fi')hp{fi')dfi', (A.7) 

Jb 



and then we verify that 

TpDpifi, n') = T^Dpifi, fi'). 



(A.8) 



The first result (IA.7l i is easily verified by integration by parts: 
for any two functions <^{fi) and r]{ii), it may be shown that, 
whether b = cos Q or b — —1, 



[b )(i ):rT' 

L dfi dfi 

+ Lp{Lp + 3)fi^ Cv 



iTpC)vdfi. 



(A.9) 



To verify the second result ( IA.8> we use the Laplace- 
Beltrami identity V,^jP;,„ = -l(l + l)Pi^ dPahlen & Tromd 
tl998) to write 



{rp~r;)Dp{fi,fi') 



(A. 10) 



„ ^'2) Cl,nPlm{f^)Plm{f^') 

I— nip 

X [l{l + l)-Lp{Lp + 3)] 

- 2fi{i - 1?) Chn-^PlMPlm{^i') 



+ 2/x'(l - /i") ClmPlM-^Plmifi'). 

I— rap 

Two well-known Legendre identities help us simplify 
the above; a derivative identity and a recursion relation 
(Dahlen & Tromp 19 98): 



dPh 



dfi 

fJ.Plni 



{I + l)fiPlm -{l-m+ \)Pl+lm 

(/ - TO + l)P; + i,n + + m)Pi^i,n 
21 + 1 



(A. 11a) 
(A. lib) 



Applying eq. ( IA.lla> . then eq. ( lA.l lb> . transforms 
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eq. ( IA.10> into 

{%-r;)D^{^i,ii') 



(A. 12) 



I— nip 

X [(I - 2){l + 1) - LpiLp + 3)] 

Lp 

+ Cl^[Pl+2MPlm{^i') - PlMPl+2m{f^')] 



l—rn„ 



X 2 



{l~m + l){l-m + 2) 



21 + 3 

The Legendre derivative identity of eq. ( lA.l lat can be ma- 
nipulated to yield a formula of the Christoffel-Darboux type 
(kSimons et al.,2006i) : 



(A.13) 



l—rrLr, 



[PL^+2 7nin)PLpm{fJ-') - -Pl^ m (Ai)^ip+2 m (m')] 

^ (Lp-m + 2)! 
{2Lp + 3){Lp + m)\' 



Inserting eq. ( IA.13t into eq. ( IA.12t yields: 



(A. 14) 



(/i' - m") ClmPlMPlmili') 

I— nip 

X [{l^2){l + l)-Lp{Lp + 3)\ 

+ (m^ - 2(21 + 1) Y' C„™P„™(M)^m„(/i')- 

I — nip n—nip 

Interchanging the order of summation and relabeling the 

sums. 



iTp~-T;)Dpifi,fi') = 

Lp 

- m") Y' ClraPlMPlmifJ-') 



(A.15) 



{I - 2){l + 1) - Lp{Lp + 3) + 2(2n + 1) 



n—l 



The term in square brackets always vanishes; in other words, 
{Tp -_72)Dp{^, /./) ~ and the commutation relation of 
eq. iA.bt is confirmed. Since commuting operators have the 
same eigenfunctions, we can find the spacelimited, fixed- 
order eigenfunctions h{d) or by solving the integral 
equation eq. ( IA.2L or, equivalently, by solving the differen- 
tial equation 



Tphpiii) = x/ip(/x), 



(A. 16) 



on the domain of the double polar cap, where x 7^ A is the 
associated eigenvalue. The operator Tp is Sturm-Liouville, 
i.e., it is of the form 



{ph'Y -qh + xph = 0, cos e < ^ < 1, 



(A. 17) 



,q(A*)=m2(l-^2)-l(^2 



where p(^) = (/x^ — — /i^) 

b^) — Lp{Lp + 3)/i, = 1, and the prime denotes dif- 
ferentiation with respect to /i. Since Tp is a Stu rm-Liouville 
operator, the eigenvalue spectrum of eq. ( IA.16> is simple: the 
spacing between adjacent fixed-order eigenvalues is roughly 
equant, as we illustrated in Figure|9] 



APPENDIX B: THE GRUNBAUM MATRIX 



The domains of eqs iA.2\ or (IA.16> . originally restricted to 
the area contained within the double polar caps, may be ex- 
tended to the entire sphere, < cos 6 < TT,hy writing 



-1 < ^ < 1. 



(B.l) 



The unknown functions, see eq. 
limited rather than spacelimited: 



9p 



Y 9lmXlm{0), 



h8}, must now be band- 



(B.2) 



where, as in eq. ilOi . depending on the requested order m 
and concentration bandwidth L of the solutions, irip is ei- 
ther m or m + 1, Lp is either L or L — 1, and the primed 
summation skips every second integer. As a result, gim 
requires no further identification. Inserting the representa- 
tion of eq. iB.2\ into eq. iB.l\ . multiplying both sides by 
27rsin6'X;„i(6'), integrating over all colatitudes < 9 < tt, 
and invoking the orthogonality eq. Q easily transforms 
eq. iB.H into an algebraic eigenvalue equation for the co- 
efficients of the unknown functions g: 

TpQp = XQp, (B-3) 
where we define 



TP, =2tt f Xi^{TpXi,„,)?<ine 



dO. 



(B.4) 



To avoid clutter, we have changed the subscript p indicating 
the parity of the solutions into a superscripted p on the ma- 
trix elements Tfi, . These are indexed by the integer degrees 
I and Since in Tp the only degrees involved range from 
nip to Lp, with every second degree skipped, when m = 0, 
the first element of the matrix Tg is thus Tqq, the second TIq, 
and so on, whereas the first and second elements of To are 
T°i and T3Q, and so on, respectively. The matrix whose di- 
agonalization leads to the even functions is thus given by 



t: 



\ 



V r; 



T\ 



T: 



(B.5) 



Its dimensions are [(Le — "i)/2 + 1] x [{L^ — m)/2 + 1]. The 
[{Lo - m - l)/2 + 1] X [(Lo - m - l)/2 + l]-dimensional 
matrix yielding the odd functions go is 



/ T° 

' rn 



m+lm+l 

??l+3 ?7Trf-l 



7T1-I-3 7TH-3 



T 



mr\-3 Lo 



^LoLo 



(B.6) 
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The dimensions of the combined matrix 
T' = diag (Te,To) 



(B.7) 



are thus [L — m + 1) x {L — m + 1), as expected. Such a 
matrix can be thought of as a permutation of a "full form" 



T 



Tll 



(B.8) 



where the degrees are arranged in the correct order, without 
skipping entries. It is this matrix T that commutes with the 
matrix D of eq. i56\ . which is also in the form 



D 



Dlt. 



Dr, 



(B.9) 



D 



LL 



(B.IO) 



and whose elements are given by eq. (I56> . 

A;' = 27r / XiraXiiryiSmedO. 

JO 

Deriving the form of the matrix entries Tf^, o f eq. re- 
quires the Griinbaum operator 7p of eq. ( lA.U as a function 
of colatitude, i.e. 



CDS' 



ml 



2 cos 6' sin 6* — 
au 



— Lp{Lp + 3) cos^ 



(B.ll) 



Evaluating eq. iB.4\ is perhaps not as pedestrian as by "sim- 
ply reading off directly the inner products" as claimed by 
[Griinbaum et al. ( 1982), but since only the result is of any 
practical consequence here, we will simply state it: 



1 
3 



1) cos^ e 

-2)^ + 1) 



21 + 3 
- Lp{L,p 

1(1 + 1) 



[{1 + If - m^] 
+ 3)] 



rpP 



(2l + 3){2l - 1) 
[lil + 'i)-LpiLp + 3)] 
21 + 3 



(B.12a) 



(B.12b) 



[(/ + 2)^ 



{2l + 5){2l + l) 



= otherwise. 



(B.12c) 

This is the result given as eq. ( I86> in the main text. The ele- 
ments again are specified in terms of the degree and not by 
the matrix index. Both T,. and To are thus not only real and 
symmetric, but also tridiagonal. The important result is that 
the coefficients of the even or odd optimally concentrated an- 
tipodal polar cap eigenfunctions ge{S) or go{S) both only re- 
quire the numerical diagonalization of a tridiagonal matrix. 
Both matrices have analytically prescribed elements, and a 
spectrum of eigenvalues that is guaranteed to be simple. 



